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Abstract 

Many materials such as martensitic or ferromagnetic crystals are observed 
to be in metastable states exhibiting a fine-scale, structured spatial oscillation 
called microstructure; and hysteresis is observed as the temperature, bound- 
ary forces, or external magnetic field changes. We have developed a numerical 
analysis of microstructure and used this theory to construct numerical meth- 
ods that have been used to compute approximations to the deformation of 
crystals with microstructure. 
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1. Introduction 

Martensitic crystals are observed to be in metastable states that can be mod- 
eled by local minima of the energy [1, 2, 11, 17, 19, 25, 33, 36] 

£{y) = / < M^2/( a; ); ®{ x )) dx + interfacial energy + loading energy, (1-1) 
Jn 

where f2 C R 3 is the reference configuration of the crystal, y(x) : Q — > R 3 is the 
deformation that may be constrained on the boundary dCl, and 9(x) : — > R is the 
temperature. The frame-indifferent elastic energy density <f>(F,ff) : M 3x3 x R — > R 
is minimized at high temperature 9 > 9t on SO(3) and at low temperature 9 < 9t 
on the martensitic variants U = SO(3)Z7i U • • • U SO(3)U N where the U l £ R 3x3 are 
symmetry-related transformation strains satisfying 

{RfU 1 R i :R i eG} = {U 1 , ...,U N } 



*This work was supported in part by AFOSR F49620-98-1-0433, by NSF DMS-0074043, by ARO 
DAAG55-98- 1-0335, by the Caltech CIMMS, and by the Minnesota Supercomputer Institute. 

'School of Mathematics, University of Minnesota, Minneapolis, Minnesota 55455, USA. E-mail: 
luskin@math.umn.edu 



708 



Mitchell Luskin 



for the symmetry group Q of the high temperature (austenitic) phase. The loading 
energy above results from applied boundary forces. 

Microstructure occurs when the deformation gradient oscillates in space among 
the SO(3)£/j to enable the deformation to attain a lower energy than could be 
attained by a more homogeneous state [1,25]. The simplest microstructure is a 
laminate in which the deformation gradient oscillates between RJJi £ SO(3)Ui 
and RjUj € SO(3)[/j for i ^ j in parallel layers of fine scale, but more complex 
microstructure is observed in nature and is predicted by the theory [2, 25]. 

We have developed numerical methods for the computation of microstruc- 
ture in martensitic and ferromagnetic crystals and validated these methods by 
the development of a numerical analysis of microstructure [4,6,12,14,16,22,25- 
28]. Related results are given in [9,10,15,21,24,31,32,34]. For martensitic crys- 
tals, we have given error estimates for stable quantities such as nonlinear integrals 
J Q f(x,Vy(x)) dx for smooth functions f(x,F) : Q x R 3x3 — > R and for the local 
volume fractions (Young measure) of the variants SO(3)Ui even though pointwisc 
values of the deformation gradient are not stable under mesh refinement. 

To model the evolution of metastable states, we have developed a computa- 
tional model that nucleates the first order phase change since otherwise the crystal 
would remain stuck in local minima of the energy as the temperature or bound- 
ary forces are varied [8]. Our finite element model for the quasi-static evolution of 
the martensitic phase transformation in a thin film nucleates regions of the high 
temperature phase during heating and regions of the low temperature phase during 
cooling. 

Graphical images for the computations of microstructure and phase transfor- 



mation described in this paper can be found at http://www.math.umn.edu/~luskin 



and in the cited references. A more extensive description of microstructure and its 
computation can be found at the above website as well as in the selected references 
at the end of this paper. 



2. Numerical analysis of microstructure 

Martensitic crystals typically exist in metastable states for time-scales of tech- 
nological interest. Many important analytic results have been obtained for math- 
ematical models for martensitic crystals, especially for energy-minimizing defor- 
mations with microstructure [1,2,11,20,30,35]. These results and concepts for 
energy-minimizing deformations should also have a role in the analysis of metasta- 
bility [18,29]. Similarly, we have developed a numerical analysis of microstructure 
[4, 6, 12, 14, 16, 22, 25-28] for which results have been obtained primarily for the ap- 
proximation of energy-minimizing deformations that we think also give insight and 
some validation for the investigation of metastability by computational methods. 

We give here a summary of the numerical analysis of martensitic microstruc- 
ture that we have developed for temperatures 9 < 9t for which the energy density 
4>(F, 9) is minimized on the martensitic variants U — SO(3)C/i U • • ■ U SO(3)L0v- 

We assume that the energy density (f>(F, 9) is continuous and satisfies near the 
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minimizing deformation gradients hi the quadratic growth condition given by 

(f>{F,9) > n\\F -ir{F)\\ 2 for all F G M 3x3 , (2.1) 
where fj, > is a constant and it : R 3x3 — ► U is a projection satisfying 



|F-7r(F)|| = min \\F - Gil for all Fe 



p3x3 



We also assume that the energy density 4>(F, 8) satisfies the growth condition for 
large F given by 



<j)(F, 0) > C^FY - Co for all F £ 



p3x3 



where Co and C\ are positive constants independent of F £ R 3x3 where p > 3 to 
ensure that deformations with finite energy are uniformly continuous. 
We can then denote the set of deformations of finite energy by 

W* = {y e C{fl;R 3 ) : [ (/>{Vy(x), 9) dx < oo }, 



and we can define the set A of admissible deformations to be 

A = { y e : y{x) = y (x) for all x £ dfl}. (2.2) 

Since we assume that the set of admissible deformations A is constrained on the 
entire boundary dtt, we can neglect the loading energy in i|l.lfl . We will also set 
the interfacial energy to be zero in this section so as to consider the idealized model 
for which the length scale of the microstructure is infinitesimally small. For the 
theorems below, we assume boundary conditions compatible with a simple laminate 
mixing QUi for Q € SO (3) with volume fraction A and Uj with volume fraction 1 — A, 

y (x) = [XQUi + (1 - X)Uj] x for all x ett, 

where for aeK 3 and net 3 , with a, n ^ 0, we have the interface equation [1, 19, 25] 

QUi = Uj + a® n. 

We consider the finite element approximation of the variational problem 

m££(y) 

given by 

inf £{y h ) 

VhSAh 

where Ah is a finite-dimensional subspace of A defined for h £ (0, ho] for some 
ho > 0. The following approximation theorem for the energy has been proven for 
the Pk or Qk type conforming finite elements on quasi-regular meshes, in particular 
for the Pi linear elements defined on tetrahedra and the Q\ trilinear elements defined 
on rectangular parallelepipeds [4, 10, 22, 25-27]. 
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Theorem 2.1. For each h € (0, ho], there exists yh € Ah such that 

£{y h ) = min £{z h ) < Ch 1 ' 2 . (2.3) 

We next define the volume fraction that an admissible deformation y 6 A is 
in the fc-th variant SO(3)[/ fe for k £ {1, ... , TV} by 

measf2 fc (?/) 

T fc 2/) = FT - 

mcas i I 

where 

Q fe (y) HieSl: 7r(Vy(x)) e SO(3)C/ fc }. 

The following stability theory was first proven for the orthorhombic to mono- 
clinic transformation (N — 2) [26] and then for the cubic to tetragonal transfor- 
mation (N = 3) [22]. The analysis of stability is more difficult for larger N since 
the additional wells give the crystal more freedom to deform without the cost of 
additional energy. In fact, for the tetragonal to monoclinic transformation (N = 4) 
[6], the orthorhombic to triclinic transformation (N = 4) [16], and the cubic to or- 
thorhombic transformation (N = 6) [4] we have shown that there are special lattice 
constants for which the laminated microstructure is not stable. Error estimates are 
obtained by substituting the approximation result 12.3f) in the following stability 
results. 

In each case for which we have proven the approximation of the microstructure 
to be stable, we have derived the following basic stability estimate for the approxi- 
mation of a simple laminate mixing QUi and Uj which bounds the volume fraction 
that y S A is in the variants k =/= i, j 

T k {y) < C (£(y)i + £(y)) for all k?i,j and y e A. (2.4) 

For the theorems that follow, we shall assume that the lattice parameters are such 
that the estimate (|2.4|l holds. 

The following theorem gives estimates for the strong convergence of the pro- 
jection of the deformation gradient parallel to the laminates (the projection of the 
deformation gradient transverse to the laminates does not converge strongly [25]), 
the strong convergence of the deformation, and the weak convergence of the defor- 
mation gradient. 

Theorem 2.2. (1) For any Hiel 3 such that w ■ n — and \w\ — 1, we have 
the estimate for the strong convergence of the projection of the deformation gradient 

J | (Vy(x) - Vy (x)) w\ 2 dx < C (s(y) + for all y e A. 

(2) We have the estimate for the strong convergence of the deformation 

J \y(x)- y (x)\ 2 dx <c(£{y)+£(y)^) for all y e A. 
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(3) For any Lipshitz domain lu <Z Q, there exists a constant C = C{uj) > 
such that we have the estimate for the weak convergence of the deformation gradient 



/ 



{S7y(x)-S7y (x))dx 



<c(£(y)i +£(y)i) for all ye A. 



For fixed i,j with i ^ j we define a projection operator mj : R 3x3 — » SO(3)/7iU 
SO(3)Uj- by 

\\F-Tnj{F)\\ ={\\F-G\\ :GeSO(3)C/ i USO(3)f/ J } for all F e M 3x3 , 

and the operators 6 : R 3x3 -» 5*0(3) and n : R 3x3 -» {QUi,Uj} by the unique 
decomposition 

TTij(F) = 0(_F)II(.F) for all F g R 3x3 . 

The next theorem shows that the deformation gradients of energy-minimizing 
sequences must oscillate between QUi and Uj. 
Theorem 2.3. We have for all y g A that 

J \\Vy(x) - n(Vy(x))|| 2 dx < C (s(y)+E{y)i) . 

We now present an estimate for the local volume fraction that a deformation 
y g A is near QUi or Uj. To describe this, we define the sets 

^ P {y) = {xeuj: U(Vy(x)) = QU t and \\Vy{x) - QUi\\ < p}, 

<4(V) = i x e uj : n(Vy(x)) = (7, and ||Vy(x) - f/J < p}, 

for any subset w g Q, p > 0, and y g A The next theorem demonstrates that 
the deformation gradients of energy-minimizing sequences must oscillate with local 
volume fraction A near QUi and local volume fraction 1 — A near Uj. 

Theorem 2.4. For any Lipshitz domain uj C ft and for any p > 0, there 
exists a constant C — C(uj, p) > such that for all y g A 



measw*(y) 



- A 



meas uj 



meas uj 3 p (y) 



mcasw 



-(1-A) 



< C 



(£(y)*+£(y)->). 



We next give an estimate for the weak stability of nonlinear functions of de- 
formation gradients. 

Theorem 2.5. We have for all f : O x R 3x3 -^R andy eA that 



J 

Jn 



m 

where 



{f(x, Vy(a:)) - [A/(x, Qt/ 4 ) + (1 - X)f(x, Uj)}} dx < C\\f\\ v E{y)* + £{yY 



\l 



with Zf : Ct 



/ {(esssup||V f /(x,F)||) 2 + \X7z f (x)n\ 2 + z f (x) 2 } dx < oo 
Jn 

R defined by 

zj(x) — f(x, QUi) — /(x, Uj) for all ie!l. 
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3. A computational model for martensitic phase 
transformation 

We have developed a computational model for the quasi-static evolution of 
the martensitic phase transformation of a single crystal thin film [8] . Our thin film 
model [7] includes surface energy, as well as sharp phase boundaries with finite 
energy. The model also includes the nucleation of regions of the high temperature 
phase (austenite) as the film is heated through the transformation temperature and 
nucleation of regions of the low temperature phase (martensite) as the film is cooled. 
The nucleation step in our algorithm is needed since the film would otherwise not 
transform. 

For our total-variation surface energy model, the bulk energy for a film of 
thickness h > with reference configuration flh = ft x (—h/2, h/2), where ft C K 2 
is a domain with a Lipschitz continuous boundary dQ, is given by the sum of the 
surface energy and the elastic energy 

k [ |D(Vu)|+ / (j)(\7u,9)dx, (3.1) 

where J n \D(Vu)\ is the total variation of the deformation gradient [7] and k is a 
small positive constant. 

We have shown in [7] that energy-minimizing deformations u of the bulk energy 
1)3. ljl are asymptotically of the form 

u(xi,x 2 ,x 3 ) = y(x 1 ,x 2 ) + b{xi,x 2 )xz +0(3:3) for (xi,x 2 ) G ft, x 3 G h/2), 

(which is similar to that found for a diffuse interface model [3]) where (y, b) mini- 
mizes the thin film energy 

8 (y, b, 9) = K ( [ \D(Vy\b\b)\ + V2 f \b- b \) + f 0(V#, 9) dx (3.2) 
\Jn Jon / Jn 

over all deformations of finite energy such that y — yo on d£l. The map b describes 
the deformation of the cross-section relative to the film [3]. We denote by (Vy|6) G 
M 3x3 the matrix whose first two columns are given by the columns of Vy and the 
last column by b. In the above equation, j n \D(Wy\b\b)\ is the total variation of the 
vector valued function (Wy\b\b) : il — > M 3x4 . 

We describe our finite element approximation of l|3.2|) by letting the elements 
of a triangulation r of ft be denoted by K and the inter-element edges by e. We 
denote the internal edges by e C Vt and the boundary edges by e c dQ. We define 
the jump of a function if> across an internal edge e C ft shared by two elements 
Ki,K 2 G r to be 

Me = 1pe,Ki ~ 1pe,K 2 , 

where tpe,Ki denotes the trace on e of ^P\ki, and we define ?/>| e to be the trace on 
e for a boundary edge e C dVt. Next, we denote by Vi(t) the space of continuous, 
piecewise linear functions on ft which are linear on each K G r and by Vq(t) the 
space of piecewise constant functions on ft which are constant on each K G t. 
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Finally, for deformations (y, b) <G V\{t) x Voir) and temperature fields 9 € Vq(t) 
the energy l|3.2|) is well-defined and we have that 



/ \D(Vy\b\b)\ + V2 f \b-b \ 
Jn Jon 



+ / <f>(Vy\b,6)dx 



<> 



= « £|[(V#|h)] e ||e| + V2 £ \ b U-bo\e \e\\+J2<t>((Vy\b,e)\ K ) \K\, 

where | • | denotes the euclidean vector norm, \e\ denotes the length of the edge e, 
\K\ is the area of the element K, and 



(V#|6)] e | = (|[Vy] e | 2 + 2|[6] £ 



1/2 



The above term is not differentiable everywhere, so we have regularized it in our 
numerical simulations. 

Since martensitic alloys are known to transform on a fast time scale, we 
model the transformation of the film from martensite to austcnite during heat- 
ing by assuming that the film reaches an clastic equilibrium on a faster time 
scale than the evolution of the temperature, so the temperature 9(x,t) can be 
obtained from a time-dependent model for thermal evolution [8]. To compute the 
evolution of the deformation, we partition the time interval [0, T] for T > by 
= to < t\ < ■ ■ ■ < <tL=T and then obtain the solution (y(te),b(ti)) e A T 
for I = 0, . . . , L by computing a local minimum for the energy £(v,c,6(te)) with 
respect to the space of approximate admissible deformations 

At = {(v,c) e V x {t) x V {r) : v = y on Oft}. (3.3) 

Since the martensitic transformation strains IA cR 3x3 are local minimizers of the 
energy density <j>(F, 9) for all 9 near 9t , a deformation that is in the martensitic 
phase will continue to be a local minimum for the bulk energy £(v,c,9(t)) for 
9 > 9t- Hence, our computational model will not simulate a transforming film 
if we compute {y(t/),b(tt)) £ A T by using an energy-decreasing algorithm with 
the initial state for the iteration at tg given by the deformation at tt—\, that is, 
if {y^'iti), &l°l(tf)) = (y(tt-i),b(tt-i)). We have thus developed and utilized an 
algorithm to nucleate regions of austenite into (y(tt-i),b(tt-i)) € A T to obtain an 
initial iterate (y^(te),b^(te)) S A T for the computation of (y(ti),b(tg)) S A T . 

We used an "equilibrium distribution" function, P(9), to determine the prob- 
ability for which the crystal will be in the austenitic phase at temperature 9 and we 
assume that an equilibrium distribution has been reached during the time between 
tl-\ and tg. The distribution function P{9) has the property that < P{9) < 1 
and 

P{9) -> as 9 -> -oo and P{9) -> 1 as 6 oo. 

At each time te, we first compute a pseudo-random number cr(K, i) € (0, 1) on 
every triangle K € t, and we then compute (y^(tg), b^'fa)) € A T by (xk denotes 
the barycenter of K): 
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1. Ifa(K,e) <P(0(x K ,te)) and {Vy(x K ,te-i)\b(x K ,te-i),e(x K ,te)) is in 
austcnite, then set 

(|/ [0I (t/),& [0] (t/)) = (y(U-i), b(t e -i)) on K. 

2. Ifa(KJ) <P(6{x K ,t t )) and (V?/^, *<-i)|6(^. U-i),9(x K , U)) is in 
martensite, then transform to austenite on K. 

3. Ua(K,£) >P(6{x K ,t e )) and {Vy(x K , t t -{)\b(xK, te-i),6{x K , t t )) is in 
austenite, then transform to martensite on K. 

4. Ua(KJ) >P(6(x K ,te)) and (Vy(x K ,te-i)\b(x K ,te-i),0(xK,te)) is in 
martensite, then set 

(y^(te),b^(te)) = (y(t/_i),6(t,_i)) on X. 

We have shown in [8] for a thin film of a CuAINi alloy in the "tent" configura- 
tion that we can compute the nucleation above by setting y^(te) = y(t(,-i) G V\{t) 
and by updating the piecewise constant b^°\te) G Pq{t) by 



y,i(xK,k-i) x y,2(x K ,te-i) 

to nucleate austenite and 



b lvl (x K ,t e ) = —, ^ on K 

\y tl (x K ,k-i) x y !2 {x K ,te-i)\ 



,[oi/ , \ y,i(xK,te-i) x y, 2 (x K ,tt-i) 

b [ \XkM) =7 1 7 3^ \ 7 7 TT on 

x y,2( x K,te-i)\ 

to nucleate martensite. 

We then compute (y(ti),b(tej) G A T by the Polak-Ribiere conjugate gradient 
method with initial iterate (y^(tf),b^(te)) G A T . We have also experimented 
with several other versions of the above algorithm for the computation of b^(tg). 
For example, the above algorithm can be modified to utilize different probability 
functions P(9) in elements with increasing and decreasing temperature. We can also 
prohibit the transformation from austenite to martensite in an element in which 
the temperature is increasing or prohibit the transformation from martensite to 
austenite in an clement for which the temperature is decreasing. 
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